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Abstract 

We present a refinement of the Spectral Method by incorporating an optimization 
method into it and generahze it to two space dimensions. We then apply this Refined 
Spectral Method as an extremely accurate technique for finding the bound states of the 
two dimensional time-independent Schrodinger equation. We first illustrate the use of 
this method on an exactly solvable case and then use it on a case which is not so. This 
method is very simple to program, fast, extremely accurate {e.g. a relative error of 10^^^ 
is easily obtainable in two dimensions), very robust and stable. Most importantly, one 
can obtain the energies and the wave functions of as many of the bound states as desired 
with a single run of the algorithm. 

PACS numbers: 02.70.Hm, 03.65. Ge 



1 Introduction 

Eighty years after the birth of quantum mechanics [1], the Schrodinger's famous equation 

still remains a subject for numerous studies, aiming at extending its field of applications and 

at developing more efficient analytic and approximation methods for obtaining its solutions. 

There has always been a remarkable interest in studying exactly solvable Schrodinger equations 

which has been found for only a very limited number of potentials, most of them being classified 
*Email: pedram@sbu.ac.ir 
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already by Infeld and Hull [2] on the basis of the Schrodinger factorization method [3] , which 
in turn appeared to be a rediscovery of the formalism stated nearly 120 years ago by Darboux 
[4]. However, a vast majority of the problems of physical interest do not fall in the above 
category when we formulate a more or less realistic model for them. Then we have to resort 
to approximation techniques which can be analytic or numeric. The Schrodinger equation can 
always be solved numerically, which nowadays seems elementary, in view of the immensely 
increased computational power. However, even in this simplest case, the success of applying 
any direct numerical integration method depends on the quality of initial guesses for the 
boundary conditions and energy eigenvalues. Moreover, one usually encounters difficulties with 
the intrinsic instabilities of typical problems, and rarely with the existence of actual solutions 
which posses rapid oscillation. The need for evermore accurate and efficient numerical methods 
for solving problems of physical interest have stimulated development of more sophisticated 
integration approaches, e.g. embedded exponentially-fitted Runge-Kutta [5] and dissipative 
Numerov-type [6] methods, as well as interesting techniques, such as a relaxation approach [7] 
based on the Henyey algorithm [8], an adaptive basis set using a hierarchical finite element 
method [9], and an approach based on microgenetic algorithm [10], which is a variation of 
a global optimization strategy proposed by Holland [11]. Most of these methods are either 
completely designed for the one-dimensional cases or optimized so. Few general methods 
are readily available for higher dimensional cases, e.g Finite Element Method (FEM), Finite 
Difference Method (FDM), Relaxation Method, Spectral Element Method. 

Here we extend the Refined Spectral Method (RSM), which is introduced in Ref. [12], 
as a numerical method to solve the higher dimensional schrodinger equations. There, we 
first refined the Spectral Method (SM) [13] for one-dimensional cases by incorporating an 
optimization procedure into it, and then tested the results obtained by our method against the 
corresponding values of an exactly solvable case. We showed that this method can be extremely 
accurate, {e.g. errors of order 10~^^°), and has the following advantages: It is very simple, 
fast, very robust and stable, i.e. it does not have the instability problems due to the usual 
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existence of divergent solutions of most physical problems. These problems usually produce 
difficulties for the spatial integration routines such as FDM and FEM. Finally, and perhaps 
most importantly, we can obtain the wave functions and energies of as many of the bound 
states as desired with a single run of the algorithm. Spectral Method, consists of first choosing 
a complete orthonormal set of eigenstates of a, preferably relevant, hermitian operator to be 
used as a suitable basis for our solution. For this numerical method we obviously can not 
choose the whole set of the complete basis, as these are usually infinite. Therefore we make 
the approximation of representing the solution by a superposition of only a finite number of 
the basis functions. By substituting this approximate solution into the differential equation, 
a matrix equation is obtained. The energies and expansion coefficients of these approximate 
solutions could be determined by the eigenvalues and eigenfunctions of this matrix, respectively. 
In the Spectral Method the concentration is on the basis functions and we expect the final 
numerical solution to be approximately independent of the actual basis used. Moreover in 
this method, the refinement of the solution is accomplished by choosing a larger set of basis 
functions, rather than choosing more grid points, as in the numerical integration methods. For 
more detailed explanation on this subject, in particular different branches of SM, including 
the commonly used Pseudo-Spectral Method, and its historical development see for example 
Ref. [13]. For an interesting appfication of this method to the double well potential see for 
example Ref. [14]. 

The remainder of this paper is organized as follows. In Section 2, wc present the underlying 
theoretical bases for the formulation of the RSM and introduce our optimization procedure in 
two-dimensions. In Section 3, we first use this method for the 2D Simple Harmonic Oscillator 
(2D-SH0), which is an exactly solvable problem, to illustrate and test the method. In Section 
4, we apply this method to an interesting 2D problem which could be relevant to QCD and is 
not exactly solvable. In Section 5, we state our conclusions. 
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2 The Refined Spectral Method 

Let us consider the 2-D time-independent Schrodinger equation, 

where m, U{x,y), and E stand for the reduced mass, potential energy, and energy, respec- 
tively. Obviously, This is an eigenvalue problem with eigenfunction ilj{x,y) and eigenvalue 
E. Throughout this paper, we only examine the bound states of this problem, i.e. the states 
which are the square integrable. Therefore the general eigenvalue problem that we want to 
solve can be cast in the form of a linear elliptic PDE one that can be written as, 

((P^j{x,y) d^tl){x,y)\ , ?/ x , / x ,/ x 



where. 



dx"^ dy'^ 



f{x,y) = '^U{x,y), 



The configuration space for most physical problems are defined by — oo < x,y < oo. We make 
the approximation of constraining the domain to —Lx/2 < x < Lx/2 and —Ly/2 < y < Ly/2. 
As shall be explained later, first of all, this constraining of the domain is absolutely crucial 
for our method, and secondly does not necessarily pose a loss of accuracy: The use of a 
finite domain is necessary since we need to choose a finite subspace of a countably infinite 
basis. Moreover, since the bound states have compact support, a finite region suffices, and 
the choices of and Ly are in fact the essential part of our optimization procedure. As 
mentioned before, any complete orthonormal set can be used for the SM. We use the Fourier 
series basis as an example. For this particular basis, we find it convenient to shift the domain 
to < X < La; and < y < Ly. In particular, we need to shift the potential energy functions 
also. This means that we can expand the solution as, 

,/ \ A ■ f'm-Kx\ . fmry\ 

^{X: y)= 2^ Am,n Sm I — — Sm -— . (4) 



We can also make the following expansion, 

f{x, y)il^{x, y) = Y. Brn,n sin ( ) sin I ^ ) , (5) 

where Bm,n are coefficients that can be determined once f{x, y) is specified. By substituting 
Eqs. (4,5) into Eq. (2) and using the differential equation of the Fourier basis we obtain, 
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E 



. / rrnix \ . j nny 



sm 



Because of the linear independence of sin and sin , every term in the summation 

must satisfy, 

/ TJITT \ ^ / 77/ TT \ i 

J V^^^J j ~ ^ ^'''^ 

It only remains to determine the matrix B. Using Eq. (5) and Eq. (4) we have, 

^ /mnx\ ( nny\ . i:, . / rmix\ fmiy\ , . 

Bm,n sm — — sm -— = 2^ Am,nf{x, y) sm — — sm -— . (8) 

m,n \ J^x y \ J m,n V i^a; / \ ^y J 

By multiplying both sides of the above equation by sin sin and integrating over 

the x, |/-space and using the orthonormality condition of the basis functions, one finds, 

Bm,n ^ ^ Cm,m' ,n,n' -^m' ,n' ) (9) 
m',n' 

where, 

4 /"^^ /"^j/ . fmTix\ . (miy\ ?, , . (m''Kx\ . (n'Tiy\ , , 
<^m,m',n,n' = (i-^) / / sm — — sm -— f[x,y) sm — — sm — — dxdy.[lQ) 
L^Ly Jo Jo \ J \ J \ J \ ^y J 

Therefore we can rewrite Eq. (7) as. 



m' ,n' 



^ J \^y , 

It is obvious that the presence of the operator f{x,y) in Eq. (2), leads to nonzero coefficients 
Cm,m',n,n' in Eq. (11), which in principle could couple all of the matrix elements of ^4. It is easy 
to see that the more basis functions we include, the closer our solution will be to the exact 



one. We select a finite subset of the basis functions i. e. the first N'^ ones, by letting the index 
m and n run from 1 to in the summations. Then we replace the square matrix A with a 
column vector A' with N'^ elements, so that any element of A corresponds to one element of 
A'. With this replacement, Eq. (11) can be written as, 

DA'^sA', (12) 

where is a square matrix with {N'^) x (A^^) elements. Its elements can be obtained from 
Eq. (11). The eigenvalues and eigenfunctions of the Schrodinger equation are approximately 
equal to the corresponding quantities of the matrix D. That is the solution to this matrix 
equation simultaneously yields N"^ sought after eigenstates and eigenvalues. The only problem 
which remains is to solve the eigenvalue problem Eq. (12), and to control the round-off errors. 
This is often a serious issue for the usual spatial integration method using double precision. 
However, we can easily overcome this problem and obtain a very high precision. Using RSM 
in ID accuracies of order 100 significant digits are very easily accomplishable while in 2D 10 
significant digits are obtained using the same computation time. This can be implemented, 
for instance with MATHEMATICA, using the instruction 'Set [Precision [...,20]', for example, 
to set a precision of 20 digits for the numbers. This method, in principle, allows us to obtain 
the eigenvalues and eigenvectors with a maximum precision of 20 digits (using enough basis 
elements) . 

Now wc can introduce our optimization procedure. We are free to adjust two parameters: 
A^, the number of basis elements used and the lengths of the spatial region, and Ly. These 
lengths should be preferably larger than spatial spreading of all the sought after wave functions. 
However, if and Ly are chosen to be too large we loose overall accuracy. After fixing these 
lengths, any desired accuracy can be obtained with a suitable choice of N . As we shall show, 
the error decreases extremely rapidly as the number of basis elements is increased. However, it 
is important to note that for each A^, and Ly have to be properly adjusted. We shall denote 
these optimal quantities by Lx{N) and Ly{N). We have come up with a method to determine 
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these quantities: For a few fixed values of N we compute E{N, L^, Ly) which invariably has 
an minimum point. Therefore, all we have to do is to compute the position of these minimum 
points and compute an interpolating function for obtaining Lx{N) and Ly{N). Obviously the 
more points we choose the better our results will be. As we shall see, the addition of this 
refinement can have dramatic consequences. 

Computation of the relative error in the exactly solvable cases is straightforward. For 
example for computing the relative error of the eigenvalue, denoted by 5^, we only need to 
find the absolute value of the difference between the result and the exact one and divide by 
the latter. For cases which are not exactly solvable, we compute the difference between the 
eigenvalues for a given N and those obtained with 1, both lying on the Lx{N) and Ly(N) 
curves. Wc shall denote the error computed by this procedure 6e- We have computed Lx{N) 
and Ly{N) for all cases, and subsequently computed the eigenf unctions, eigenvalues and their 
errors using this method, and checked their vahdity in the exactly solvable case of 2D-SH0. 
Obviously to obtain consistent results we have to keep the same precision throughout the 
calculations. 

3 2D Simple Harmonic Oscillator 

In this section, for illustrative purposes, we apply RSM to find the bound states of a 2D-SH0. 
We can then readily check the validity of our whole procedure, which includes our prescription 
for finding the optimal quantities Lx{N) and Ly{N), and the overall accuracy of our results. 
The Schrodinger equation for an isotropic 2D-SH0 is, 

~dx^ dy^ l + -mw(x + y )ilj{x ,y ) = E ip{x ,y ), (13) 

where uj is the natural frequency of the Oscillator. We first shift the variables as explained 
above, and then we convert this differential equation into the following dimensionless form by 
dividing both sides by huj/2, 

-^^^ - ^^^^ + ((^ - L./2r + {y- Lj2f) ^(x, y) = E^l^^x, y), (14) 
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where x = ^''^x\y = ^J^y' ■, and E — ^E'. This differential equation is exactly solvable and 
its eigenvalues and eigenf unctions, which are all bound states, can be easily found analytically 
and are well known, 

= (i)"'%MS^e-<«^-»^)^ (16) 

-^n„n^ = 2{n^ + ny + l), no:, % = {0,1,2, ...}, (16) 



where Hn{x) denote the Hermite polynomials. Using RSM we can calculate accurately the 
energy levels and the corresponding eigenfunctions of this Hamiltonian. Here, we choose our 
optimization procedure for the ground state which will be symmetric in x and y, therefore 
Lx{N) = Lx{N) = L{N). The computation of the errors of the wave functions are analogous 
to that of the energy. We divide the configuration space into M grid points. Then, we average 
the square of the absolute value of the difference between the exact solution and that obtained 
by the RSM on the grid points, 

,2 _ ^f!j=l\^e.act{ij)-MiJ)\' e 

In the above equation we have also shown the expression for Se, for ease of reference. Figure 
1 shows the ground state energy computed using SM for the fixed value of the = 6 as a 
function of L. Note the existence of the minimum point at the exact value of the eigenvalue. 
This point determines 1/(6). We repeat this procedure for a few other values of N. After 
plotting these values we can obtain an interpolating function L{N) (Fig. 2). The optimization 
method introduce here is equivalent to the one introduced in Ref. [12], where inflection points 
determined the quantities L{N). Table 1 shows the complete results for the first 10 eigenvalues 
and eigenvectors for — 22. Several points are note worthy here. First, note the outstanding 
accuracy of 6e ~ 10~^^ for the ground state in particular, and the general good correspondence 
between Se and Se- Also note the corresponding good accuracy for S^, reported only for the 
non- degenerate cases. We did not calculate S^ for other cases because the outcome of the al- 
gorithm in each degenerate subspace causes an unpredictable linear combination of those wave 
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Figure 1: Ground state energy for 2D-SH0 versus L for N = 6, using SM in units where 
huj = 2. The position of the minimum determines 1/(6). 

functions, which is equivalent to a whole rotation in that subspace. Hence the computation 
of becomes a little complicated. Also note that the errors associated with wave functions 
symmetric in x and y are about one order of magnitude better than the asymmetric ones in 
the degenerate subspace, because we assumed this symmetry in our optimization procedure. 
In Fig. 3 we show a semi-log plot of the error for the ground state energy, obtained using 
RSM, in terms of N, all obtained using appropriate L{N). Note that the error falls off exactly 
exponentially as a function of N, a theoretical property common to all SM [13]. The exact 
matching of our computed error with this theoretical expectation is another positive sign for 
our method. In Fig. 4 we state the MATHEMATICA program for solving this problem, to 
emphasize how short our program is. We have only left out the the computation of L{N). 
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Figure 2: The dots represent the values of L computed by the method described in the text 
for different values of N. The solid hne represents the computed interpolation function L{N). 
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3,1 


10. 


10.00000000020025511 


2.00 X 10-^1 


1.90 X 10-11 




0,4 


10. 


10.00000000630282991 


6.30 X 10-1° 


5.90 X 10-1° 




4,0 


10. 


10.00000000630282991 


6.30 X 10-1° 


5.90 X 10-1° 




2,3 


12. 


12.00000000021802137 


1.81 X 10-11 


1.73 X 10-11 




3,2 


12. 


12.00000000021802137 


1.81 X 10-11 


1.73 X 10-11 




1,4 


12. 


12.00000000630309285 


5.25 X 10-1° 


4.24 X 10-1° 




4,1 


12. 


12.00000000630309285 


5.25 X 10-1° 


4.24 X 10-1° 




0,5 


12. 


12.00000003939548075 


3.28 X 10-*^ 


3.08 X 10-9 




5,0 


12. 


12.00000003939548075 


3.28 X 10-9 


3.08 X 10-9 




L(22) 


1197 

ion 



Table 1: The results for the first 21 eigenstates (out of 484) of the 2D-SH0 in units where 
hco = 2, using RSM with = 22. That is we have used 484 basis functions and L(22) = 11|^. 
Note that we obtain all the degeneracies and all with very good accuracy. Note that we have 
good correspondence between 5e and 5e- The computation time was about 697.765 seconds 
for N ^22 (using SetPrecision[..., 20] in MATHEMATICA) on a Pentium 2.4 GHz machine 
resulting in a precision of lO-i^ for the ground state energy, and it was 11.062 seconds for 
N — 18 (using default double precision) resulting in a precision of IO-12. 



11 



NN = N; 

U[x_, y_] := x'^a + y-^a 
L = L; 
i = 1; 

Nprecision = 50; 
Do[ 

, . .mTTx, rUTT-y, rinlTrx, rnlTry. 
{cl[i] = Integrate [ Sin[ ] Sin[ ] U[x-L/2, y-L/2] Sin[ ] Sin[ ] , 

{X, 0, Lx}, {y, 0, Ly}] , i = i + l}, {m, 1, NN} , {n, 1, NN} , {ml, 1, NN} , {nl, 1, NN}] ; 
Do [mat [m, n] = (2 / L) 2 cl [ (m - 1) *NN'^2 + n] + 

KroneckerDelta[IntegerPart [ (m - 1) /NN+1], IntegerPart [ (n - 1) /NN+1]] 
KroneckerDelta[Mod[m- 1, NN] + 1, Mod[n - 1, NN] + 1] 

( (IntegerPart [ (n - 1) /NN+1] tt/ L) ^^2 + ( (Mod[ (n - 1) , NN] + 1) tt/ L) -^2) , 
{m, 1, NN'^2}, {n, 1, NN'^2}]; 
Do[mat2[m, n] = SetPrecision[mat [m, n] , Nprecision], {m, 1, NN'^2}, {n, 1, NN'^2}]; 
mats = Table [mat2 [n, m] , {m, 1, NN'^2}, {n, 1, NN'^2}]; 
Value = Eigenvalues [mats] 
Vec = Eigenvectors [mats] ; 

Figure 4: MATHEMATICA commands for computing the spectrum of the Hamiltonian H — 

pI^U (x, y), in units where Ti = 1. The value for L in hne 3 should be obtained by our 
optimization procedure as described in the text. Whenever possible we evaluate the integrals 
analytically, and replace the Integrate command by its results, to increase the precision and 
save time. This has been the case for the examples presented in this paper. 
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n 


ESM 


5e 


1 


1.10822315780256 


1.19 X 10-^° 


2 


2.37863785124994 


1.16 X 10-8 


3 


2.37863785124996 


1.16 X 10-8 


4 


3.05608156130323 


2.06 X 10"^ 


5 


3.51495134040797 


1.12 X 10-6 


6 


4.09348955687600 


9.38 X 10-*^ 


7 


4.09348955687604 


9.38 X 10-6 


8 


4.75298944936096 


9.32 X 10-5 


9 


4.98538290136962 


1.75 X 10-5 


10 


5.01127928161308 


4.59 X 10-" 


11 


5.50103621623983 


7.92 X 10-* 


12 


5.50103621623990 


7.92 X 10-4 


20 


8.07437393671447 


6.64 X 10-*^ 


25 


9.27305945794927 


3.36 X 10-8 


33 


11.4718771513251 


7.24 X 10-'' 


44 


13.8662683175987 


8.33 X 10-6 


L(42) 


1553 

inn 



Table 2: The results for the eigenvalues {E^^) using RSM with N = 42, for the first 12 states 
and some other highly exited and interesting ones (as explained in the text) for the 2D-QCD 
Hamiltonian, H-p^ + x^y"^, in units where h = 1. 

4 An example which is not exactly solvable 

The dimensionless and shifted Schrodinger equation for the example that we want to solve 
here is, 

_ (f^i^ + ^^^) + « - V2)^(y - L/2mx, y) = E^l^ix, y), (18) 

where a is a positive constant. The potential in this example is sometimes called the 2D-QCD 
potential. This PDE is elliptic and not exactly solvable. Therefore, we use RSM to find its 
eigenvalues and eigenf unctions. In Table. 2 we have shown the eigenvalues for the first 12 
states and some other highly exited ones {n — {20,25,33,44}). The latter were chosen for 
their unusually high accuracy due to their symmetric form, and the fact that we have chosen 
Lx{N) = Ly{N) = L(N). Figure 5 shows the ground state wave function. Note the slight 
over extension of the wave function in the x and y directions due to the particular form of the 
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Figure 5: The wave function of the ground state of the 2D-QCD potential using N = 42 and 
L(42) = in units where h = 1. 

potential. In Fig. 6 we show the wave functions for the second, forth, fifth, and forty forth 
eigenstates. The third eigenstate is not shown because it is degenerate with the second and 
can be obtained from it by 90 degree rotation. We have decided to show the forty forth state 
because we found it interesting and its highly symmetric form produces an unusually small 
error, as explained above. 

5 Conclusions 

We have extended the Refined Spectral Method to two dimensions and used it as an extremely 
accurate method for obtaining the energies and wave functions of the bound states of the two 
dimensional time-independent Schrodinger equation. In this method a finite basis is used for 
approximating the solutions. The refinement of the method is accomplished by calculating an 
optimized spatial domain for a given number of basis elements, denoted by L{N). The criteria 
for this optimization is to minimize the energy for one of the eigenstates, usually chosen to 
be the ground state. Note that our refined method is not quite equivalent to the Rayleigh- 
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Figure 6: Some eigenfunctions of the 2D-QCD potential. Upper left: 2nd; Upper right 4th; 
Lower left: 5th, and lower right: 44th states, respectively. All the states are non-degenerate 
except the first one. Its degenerate wave function can be obtained by 90 degree rotation. The 
last one is shown because it is symmetrical, interesting, and has a relatively low error. Same 
parameters were used as for Fig. 5 
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Ritz variational method, in that we have combined the spectral method in which the wave 
function is expanded in an arbitrary basis with optimization of the spatial domain which is 
equivalent to adjusting the effective potential for the bound states. This effective potential for 
the case of the Fourier basis with confinement boundary condition is the actual potential plus 
the confining "walls" placed at the boundaries, which are separated by L[N). This refinement 
scheme usually improves the accuracy of SM drastically and this effect increases rapidly with 
A^". We applied this method to an exactly solvable problem and easily found an extraordinarily 
good agreement with the exact solutions (errors of order 10~^^ with only 22 basis functions). 
This method is very simple, fast, extremely accurate in most cases, very robust, stable, and 
there is no need to specify the boundary conditions on the slopes. Most importantly, one can 
obtain the energies and the wave functions of as many of the bound states as desired with a 
single run of the algorithm. The generalization of this method to higher dimensional cases is 
straight forward. 
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